
################
### Analysis ###
################

rm(list = ls())
gc()

library(foreign)
library(reshape2)
library(stringr)
library(stargazer)
library(effects)
library(dplyr)
library(data.table)
library(scales)
library(lfe)
library(msm)
library(sjPlot)
library(ggrepel)
library(lemon)
library(gridExtra)
library(binsreg)
library(DescTools)
library(ggpubr)
library(srvyr)
library(wCorr)

### Set working directory
setwd("~/YOURPATH/replication_files")

### Load data
load("data/cs_survey.RData")

### Recode variables
cs_survey[ , age_sqr := age * age ]

cs_survey[ , unemployed := ifelse(curnt_empl_stat_RESP %in% c("Unemployed, not looking for job","Unemployed, looking for job"), 1, 0)  ]
cs_survey[ , employed := ifelse(curnt_empl_stat_RESP %in% c("Employed"), 1, 0)  ]
cs_survey[ , student := ifelse(curnt_empl_stat_RESP %in% c("In education or training (including vocational training)"), 1, 0)]
cs_survey[ , retired := ifelse(curnt_empl_stat_RESP %in% c("Retired"), 1, 0)]

cs_survey[ highest_educ %in% c("Not completed primary education","Completed primary education") , highest_educ_v2 := "primary or less"]
cs_survey[ highest_educ %in% c("Completed high school") , highest_educ_v2 := "high school"]
cs_survey[ highest_educ %in% c("Completed vocational education of 2 years or more","Completed vocational training/education of less than 2 years") , highest_educ_v2 := "vocational training"]
cs_survey[ highest_educ %in% c("Completed college") , highest_educ_v2 := "college"]
cs_survey[ highest_educ %in% c("Completed a postgraduate degree (incl. MA, MBA)","Completed a doctoral degree (incl. J.D or M.D)") , highest_educ_v2 := "postgrad"]
cs_survey[ , highest_educ_v2 := factor(highest_educ_v2, levels=c("primary or less","high school","vocational training","college","postgrad"))]

cs_survey[ , married := ifelse(living_situation == "I live with my husband/wife/partner", 1, 0) ]
cs_survey[ , renter := ifelse(homeowner == 0, 1, 0)]

cs_survey[ , liquid_assets_USD := rowSums(cbind(ass_check_sav_USD, ass_bond_stock_USD), na.rm=T) ]   
cs_survey[ is.na(ass_check_sav_USD) & is.na(ass_bond_stock_USD), liquid_assets_USD := NA ]  

### Winsorization by country
cs_survey[ , liquid_assets_USD_v1 := lapply(.SD, function(x) Winsorize(x, probs = c(0.025, 0.975), na.rm=T) ), by=c("country_name"), .SDcols= "liquid_assets_USD" ]
cs_survey[ , ass_real_estate_USD_v1 := lapply(.SD, function(x) Winsorize(x, probs = c(0.025, 0.975), na.rm=T) ), by=c("country_name"), .SDcols= "ass_real_estate_USD" ]
cs_survey[ , tot_net_inc_USD_v1 := lapply(.SD, function(x) Winsorize(x, probs = c(0, 0.975), na.rm=T) ), by=c("country_name"), .SDcols= "tot_net_inc_USD" ]

### Create deciles: those with no assets are in the bottom groups
cs_survey[ liquid_assets_USD_v1 > 0 , liquid_assets_USD_v1_dec := lapply(.SD, function(x) ntile(x, 9)), by=c("country_name"), .SDcols= "liquid_assets_USD_v1" ]
cs_survey[ liquid_assets_USD_v1 == 0 , liquid_assets_USD_v1_dec := 0 ]
cs_survey[ , liquid_assets_USD_v1_dec := liquid_assets_USD_v1_dec + 1 ]

cs_survey[ ass_real_estate_USD_v1 > 0 , ass_real_estate_USD_v1_dec := lapply(.SD, function(x) ntile(x, 9)), by=c("country_name"), .SDcols= "ass_real_estate_USD_v1" ]
cs_survey[ ass_real_estate_USD_v1 == 0 , ass_real_estate_USD_v1_dec := 0 ]
cs_survey[ , ass_real_estate_USD_v1_dec := ass_real_estate_USD_v1_dec + 1 ]

cs_survey[ , tot_net_inc_USD_v1_dec_v2 := lapply(.SD, function(x) ntile(x, 10)), by=c("country_name"), .SDcols= "tot_net_inc_USD_v1" ]
cs_survey[ ,.N, tot_net_inc_USD_v1_dec_v2 ]

### Groups based on median split
cs_survey[ , tot_net_inc_USD_med := lapply(.SD, function(x) median(x, na.rm=T)), by=c("country_name"), .SDcols= "tot_net_inc_USD_v1" ]
cs_survey[ , liquid_assets_USD_med := lapply(.SD, function(x) median(x, na.rm=T)), by=c("country_name"), .SDcols= "liquid_assets_USD_v1" ]
cs_survey[ , ass_real_estate_USD_med := lapply(.SD, function(x) median(x, na.rm=T)), by=c("country_name"), .SDcols= "ass_real_estate_USD_v1" ]

cs_survey[ tot_net_inc_USD_v1 < tot_net_inc_USD_med & liquid_assets_USD_v1 < liquid_assets_USD_med, grp_med_liq := "Liq Assets low | Income low" ]
cs_survey[ tot_net_inc_USD_v1 < tot_net_inc_USD_med & liquid_assets_USD_v1 >= liquid_assets_USD_med, grp_med_liq := "Liq Assets high | Income low" ]
cs_survey[ tot_net_inc_USD_v1 >= tot_net_inc_USD_med & liquid_assets_USD_v1 < liquid_assets_USD_med, grp_med_liq := "Liq Assets low | Income high" ]
cs_survey[ tot_net_inc_USD_v1 >= tot_net_inc_USD_med & liquid_assets_USD_v1 >= liquid_assets_USD_med, grp_med_liq := "Liq Assets high | Income high" ]

cs_survey[ tot_net_inc_USD_v1 < tot_net_inc_USD_med & ass_real_estate_USD_v1 < ass_real_estate_USD_med, grp_med_illiq := "Illiq Assets low | Income low" ]
cs_survey[ tot_net_inc_USD_v1 < tot_net_inc_USD_med & ass_real_estate_USD_v1 >= ass_real_estate_USD_med, grp_med_illiq := "Illiq Assets high | Income low" ]
cs_survey[ tot_net_inc_USD_v1 >= tot_net_inc_USD_med & ass_real_estate_USD_v1 < ass_real_estate_USD_med, grp_med_illiq := "Illiq Assets low | Income high" ]
cs_survey[ tot_net_inc_USD_v1 >= tot_net_inc_USD_med & ass_real_estate_USD_v1 >= ass_real_estate_USD_med, grp_med_illiq := "Illiq Assets high | Income high" ]

### Groups based on tertile split
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(1,2,3) & liquid_assets_USD_v1_dec %in% c(1,2,3) , grp_dec_liq := "Liq Assets low | Income low" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(1,2,3) & liquid_assets_USD_v1_dec %in% c(8,9,10) , grp_dec_liq := "Liq Assets high | Income low" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(8,9,10) & liquid_assets_USD_v1_dec %in% c(1,2,3) , grp_dec_liq := "Liq Assets low | Income high" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(8,9,10) & liquid_assets_USD_v1_dec %in% c(8,9,10) , grp_dec_liq := "Liq Assets high | Income high" ]

cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(1,2,3) & ass_real_estate_USD_v1_dec %in% c(1,2,3) , grp_dec_illiq := "Illiq Assets low | Income low" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(1,2,3) & ass_real_estate_USD_v1_dec %in% c(8,9,10) , grp_dec_illiq := "Illiq Assets high | Income low" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(8,9,10) & ass_real_estate_USD_v1_dec %in% c(1,2,3) , grp_dec_illiq := "Illiq Assets low | Income high" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(8,9,10) & ass_real_estate_USD_v1_dec %in% c(8,9,10) , grp_dec_illiq := "Illiq Assets high | Income high" ]

### Normalize preferences 
cs_survey[ , gvt_help_unemp_num_norm := BBmisc::normalize(gvt_help_unemp_num, method = "range", range = c(0, 1))  ]
cs_survey[ , gvt_reduce_inc_diff_num_norm := BBmisc::normalize(gvt_reduce_inc_diff_num, method = "range", range = c(0, 1))  ]

### Preference groups
cs_survey[ gvt_reduce_inc_diff_num %in% c(-2,-1) , gvt_reduce_inc_diff_grp := "disagree" ]
cs_survey[ gvt_reduce_inc_diff_num %in% c(1,2) , gvt_reduce_inc_diff_grp := "agree" ]
cs_survey[ gvt_help_unemp_num %in% c(-2,-1) , gvt_help_unemp_grp := "disagree" ]
cs_survey[ gvt_help_unemp_num %in% c(1,2) , gvt_help_unemp_grp := "agree" ]


#########################
### Descriptive stats ###
#########################

d_svy_wgt <- cs_survey %>% as_survey_design(weights = wgt)

d_inc <- as.data.table(cbind(d_svy_wgt %>% group_by(country_name) %>% summarise(survey_mean(tot_net_inc_USD_v1, na.rm=T)),
                             d_svy_wgt %>% group_by(country_name) %>% summarise(survey_quantile(tot_net_inc_USD_v1, na.rm=T, quantiles = 0.5)),
                             d_svy_wgt %>% group_by(country_name) %>% summarise(survey_sd(tot_net_inc_USD_v1, na.rm=T))))

d_inc <- d_inc[ ,c(1,2,5,8)]
names(d_inc) <- c("country_name","mean","median","sd")

d_liq <- as.data.table(cbind(d_svy_wgt %>% group_by(country_name) %>% summarise(survey_mean(liquid_assets_USD_v1, na.rm=T)),
                             d_svy_wgt %>% group_by(country_name) %>% summarise(survey_quantile(liquid_assets_USD_v1, na.rm=T, quantiles = 0.5)),
                             d_svy_wgt %>% group_by(country_name) %>% summarise(survey_sd(liquid_assets_USD_v1, na.rm=T))))

d_liq <- d_liq[ ,c(1,2,5,8)]
names(d_liq) <- c("country_name","mean","median","sd")

d_illiq <- as.data.table(cbind(d_svy_wgt %>% group_by(country_name) %>% summarise(survey_mean(ass_real_estate_USD_v1, na.rm=T)),
                               d_svy_wgt %>% group_by(country_name) %>% summarise(survey_quantile(ass_real_estate_USD_v1, na.rm=T, quantiles = 0.5)),
                               d_svy_wgt %>% group_by(country_name) %>% summarise(survey_sd(ass_real_estate_USD_v1, na.rm=T))))

d_illiq <- d_illiq[ ,c(1,2,5,8)]
names(d_illiq) <- c("country_name","mean","median","sd")


d_inc[ , country_name := as.character(country_name) ]
d_liq[ , country_name := as.character(country_name) ]
d_illiq[ , country_name := as.character(country_name) ]

setkey(d_inc, country_name)
setkey(d_liq, country_name)
setkey(d_illiq, country_name)

summary_tab <- rbind(d_inc, d_liq, d_illiq)
summary_tab[ , 2:4 := lapply( .SD, function (x) round(x,0) ), .SDcols=c(2:4) ]
summary_tab[ , 2:4 := lapply( .SD, function (x) trimws(format(x, nsmall=0, big.mark=","),"both") ), .SDcols=c(2:4) ]

### Table 1 
print(xtable::xtable(summary_tab, digits=0), include.rownames=F)


### Preferences 
d_svy_wgt1 <- cs_survey[ !is.na(gvt_reduce_inc_diff_num) & !is.na(gvt_help_unemp_num)] %>% as_survey_design(weights = wgt)

share_gvt_redist <- as.data.table(d_svy_wgt1 %>% group_by(country_name, gvt_reduce_inc_diff_grp) %>% summarise(prop = survey_mean(), total = survey_total()))
share_gvt_redist <- share_gvt_redist[ !is.na(gvt_reduce_inc_diff_grp), .(country_name, gvt_reduce_inc_diff_grp, prop)]
setnames(share_gvt_redist, "gvt_reduce_inc_diff_grp", "type")

share_gvt_UI <- as.data.table(d_svy_wgt1 %>% group_by(country_name, gvt_help_unemp_grp) %>% summarise(prop = survey_mean(), total = survey_total()))
share_gvt_UI <- share_gvt_UI[ !is.na(gvt_help_unemp_grp), .(country_name, gvt_help_unemp_grp, prop)]
setnames(share_gvt_UI, "gvt_help_unemp_grp", "type")

share_gvt_comb <- cbind(share_gvt_redist[ type == "disagree",.(country_name,prop)],
                        share_gvt_redist[ type == "agree",.(prop)],
                        share_gvt_UI[ type == "disagree",.(prop)],
                        share_gvt_UI[ type == "agree",.(prop)])

names(share_gvt_comb) <- c("country","gvt_redist_disagree","gvt_redist_agree","gvt_unemp_disagree","gvt_unemp_agree")

share_gvt_comb[ , 2:5 := lapply(.SD, function(x) x*100 ), .SDcols= 2:5 ]

share_gvt_comb[ , country := as.character(country)]
setkey(share_gvt_comb, country)

### Table 2
print(xtable::xtable(share_gvt_comb, digits=2), include.rownames=F)


### Figure 3
d_mean_unemp1 <- as.data.table(d_svy_wgt1 %>% group_by(grp_med_liq) %>% summarise(survey_mean(gvt_help_unemp_num_norm, na.rm=T)))
d_mean_unemp2 <- as.data.table(d_svy_wgt1 %>% group_by(grp_med_illiq) %>% summarise(survey_mean(gvt_help_unemp_num_norm, na.rm=T)))

d_mean_red1 <- as.data.table(d_svy_wgt1 %>% group_by(grp_med_liq) %>% summarise(survey_mean(gvt_reduce_inc_diff_num_norm, na.rm=T)))
d_mean_red2 <- as.data.table(d_svy_wgt1 %>% group_by(grp_med_illiq) %>% summarise(survey_mean(gvt_reduce_inc_diff_num_norm, na.rm=T)))

d_mean_unemp1[ , dv := "Unemployment insurance"]
d_mean_unemp2[ , dv := "Unemployment insurance"]
d_mean_red1[ , dv := "Redistribution"]
d_mean_red2[ , dv := "Redistribution"]

setnames(d_mean_unemp1, "grp_med_liq", "type")
setnames(d_mean_unemp2, "grp_med_illiq", "type")
setnames(d_mean_red1, "grp_med_liq", "type")
setnames(d_mean_red2, "grp_med_illiq", "type")

d_mean_unemp1[ , asset_type := "Liquid"]
d_mean_unemp2[ , asset_type := "Illiquid"]
d_mean_red1[ , asset_type := "Liquid"]
d_mean_red2[ , asset_type := "Illiquid"]

tab_pref <- rbind(d_mean_unemp1,d_mean_unemp2,d_mean_red1,d_mean_red2)

tab_pref[ , type := gsub("Liq ","",type) ]
tab_pref[ , type := gsub("Illiq ","",type) ]

tab_pref[ type == "Assets low | Income low" , value2 := "Economically precarious\n(low assets; low income)"]
tab_pref[ type == "Assets high | Income low" , value2 := "Asset-buffered\n(high assets; low income)"]
tab_pref[ type == "Assets low | Income high" , value2 := "Income-buffered\n(low assets; high income)"]
tab_pref[ type == "Assets high | Income high" , value2 := "Truly wealthy\n(high assets; high income)"]

tab_pref[ , value2 := factor(value2, levels=c("Economically precarious\n(low assets; low income)",
                                              "Asset-buffered\n(high assets; low income)",
                                              "Income-buffered\n(low assets; high income)",
                                              "Truly wealthy\n(high assets; high income)"))]

tab_pref[ , dv := factor(dv, levels=c("Unemployment insurance","Redistribution"))]

tab_pref[ asset_type == "Liquid" & dv == "Unemployment insurance", x2 := 0.5 ]
tab_pref[ asset_type == "Liquid" & dv == "Redistribution", x2 := 0.6 ]

tab_pref[ asset_type == "Illiquid" & dv == "Unemployment insurance" , x2 := 1.2 ]
tab_pref[ asset_type == "Illiquid" & dv == "Redistribution" , x2 := 1.3 ]

tab_pref[ , lower := coef - 1.96*`_se` ]
tab_pref[ , upper := coef + 1.96*`_se` ]

cairo_pdf("output/fig_3.pdf", width=7.25, height=3.5)
(plot <- (ggplot(data=tab_pref[ !is.na(type) ], aes(y=coef, x=x2, colour=dv, fill=dv, shape=dv, label=round(coef,2)))
          + geom_errorbar(aes(ymin=lower, ymax=upper), width=0)
          + geom_point(size=2.5)
          + facet_grid(~ value2)
          + theme_bw()
          + ylab("Average policy support")
          + xlab("Asset type")
          + theme(legend.position="bottom", panel.grid.minor = element_blank(), 
                  panel.border = element_rect(colour="gray"), strip.background = element_rect(colour="gray"))
          + guides(colour = guide_legend(ncol = 2))
          + scale_color_manual(name="Social policy dimension", values=c("black","darkgreen"))
          + scale_fill_manual(name="Social policy dimension", values=c("black","darkgreen"))
          + scale_shape_manual(name="Social policy dimension", values=c(21,22))
          + scale_x_continuous(breaks=c(0.55,1.25), labels = c("Liquid","Illiquid"), limits=c(0.2,1.6))
          + coord_cartesian(ylim=c(0.63,0.77))
))
dev.off()



### -----------
### Heatmap (Appendix Figure B1)
### -----------

cty_list <- sort(as.character(unique(cs_survey$country_name)))

heat_list <- list()

for(i in 1:length(cty_list)){
  
  d <- na.omit(cs_survey[ country_name %in% cty_list[i],.(tot_net_inc_USD_v1_dec_v2,liquid_assets_USD_v1_dec,ass_real_estate_USD_v1_dec,wgt)])
  
  d_wgt <- d %>% as_survey_design(weights = wgt)
  
  d_liq <- as.data.table(d_wgt %>%
                           group_by(tot_net_inc_USD_v1_dec_v2,liquid_assets_USD_v1_dec) %>%
                           survey_tally())
  
  d_illiq <- as.data.table(d_wgt %>%
                             group_by(tot_net_inc_USD_v1_dec_v2,ass_real_estate_USD_v1_dec) %>%
                             survey_tally())
  
  d_liq[ , type := "liquid_assets"]
  d_illiq[ , type := "illiquid_assets"]
  
  setnames(d_liq,"liquid_assets_USD_v1_dec","asset_decile")
  setnames(d_illiq,"ass_real_estate_USD_v1_dec","asset_decile")
  
  d_comb <- rbind(d_liq, d_illiq)
  
  d_comb[ , cty := cty_list[i] ]
  d_comb[ , n_tot := lapply(.SD, function(x) sum(x) ), by=c("type"), .SDcols= "n" ]
  d_comb[ , share := n / n_tot ]
  
  heat_list[[i]] <- d_comb
  print(i)
}

heat_list_all <- rbindlist(heat_list)

cairo_pdf("output/si_fig_b1a.pdf", width=7, height=5.5)
(ggplot(heat_list_all[ type == "liquid_assets" ], aes(y = tot_net_inc_USD_v1_dec_v2, x = asset_decile))
  + geom_raster(aes(fill=share)) 
  + facet_wrap(~ cty)
  + scale_fill_gradient(low="grey90", high="red", name="Population\nshare") 
  + ylab("Income decile")
  + xlab("Asset decile")
  + scale_y_continuous(breaks=seq(2,10,2))
  + scale_x_continuous(breaks=seq(2,10,2))
  + theme_bw()
  + theme(legend.position="right", panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
          panel.border = element_rect(colour="gray"))
  + ggtitle("Income and liquid assets")
  
)
dev.off()


cairo_pdf("output/si_fig_b1b.pdf", width=7, height=5.5)
(ggplot(heat_list_all[ type == "illiquid_assets" ], aes(y = tot_net_inc_USD_v1_dec_v2, x = asset_decile))
  + geom_raster(aes(fill=share)) 
  + facet_wrap(~ cty)
  + scale_fill_gradient(low="grey90", high="red", name="Population\nshare") 
  + ylab("Income decile")
  + xlab("Asset decile")
  + scale_y_continuous(breaks=seq(2,10,2))
  + scale_x_continuous(breaks=seq(2,10,2))
  + theme_bw()
  + theme(legend.position="right", panel.grid.minor = element_blank(), panel.grid.major = element_blank(), 
          panel.border = element_rect(colour="gray"))
  + ggtitle("Income and illiquid assets")
)
dev.off()




##############################
### Main regression models ###
############################## 

mod1.1 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.2 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + married
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.3 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.4 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

### Appendix Table D.1
lab_covar <- c("Asset-buffered","Income-buffered","Economically-precarious",
               "Asset-buffered","Income-buffered","Economically-precarious",
               "Age","Age square","Male","Education: high school","Education: vocational training",
               "Education: college", "Education: postgrad",
               "Unemployed","Student","Retired","Employed","Homeowner","Number of children","Married")

stargazer(mod1.1,mod1.3,mod1.2,mod1.4,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3,  
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar,
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1$fitted.values),3), 
                             round(mean(mod1.3$fitted.values),3), 
                             round(mean(mod1.2$fitted.values),3), 
                             round(mean(mod1.4$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$")))

### Figure 4
reg_main_out <- data.table(rbind(cbind(mod1.1$coefficients[1:3], confint(mod1.1, level=0.95)[1:3,], confint(mod1.1, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Liquid assets", grp=rownames(mod1.1$coefficients)[1:3]),
                                 cbind(mod1.2$coefficients[1:3], confint(mod1.2, level=0.95)[1:3,], confint(mod1.2, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Liquid assets", grp=rownames(mod1.2$coefficients)[1:3]),
                                 cbind(mod1.3$coefficients[1:3], confint(mod1.3, level=0.95)[1:3,], confint(mod1.3, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.3$coefficients)[1:3]),
                                 cbind(mod1.4$coefficients[1:3], confint(mod1.4, level=0.95)[1:3,], confint(mod1.4, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.4$coefficients)[1:3])))

reg_main_out[ , 1:5 := lapply( .SD, function (x) as.numeric(as.character(x)) ), .SDcols=c(1:5) ]

reg_main_out[ , grp := gsub("grp_med_liqLiq ","",grp)]
reg_main_out[ , grp := gsub("grp_med_liq2Illiq ","",grp)]

names(reg_main_out) <- c("coef","lower95","upper95","lower90","upper90","dv","asset_type","grp")

reg_main_out[ , asset_type := factor(asset_type, levels=c("Liquid assets","Illiquid assets")) ]

reg_main_out[ dv == "gvt_help_unemp_num_norm" , dv2 := "Support for unemployment insurance" ]
reg_main_out[ dv == "gvt_reduce_inc_diff_num_norm" , dv2 := "Support for redistribution" ]
reg_main_out[ , dv2 := factor(dv2, levels=c("Support for unemployment insurance","Support for redistribution")) ]

reg_main_out[ asset_type == "Liquid assets" & grp == "Assets low | Income high" , x2 := 1 ]
reg_main_out[ asset_type == "Liquid assets" & grp == "Assets high | Income low" , x2 := 1.5 ]
reg_main_out[ asset_type == "Liquid assets" & grp == "Assets low | Income low" , x2 := 2 ]

reg_main_out[ asset_type == "Illiquid assets" & grp == "Assets low | Income high" , x2 := 1.1 ]
reg_main_out[ asset_type == "Illiquid assets" & grp == "Assets high | Income low" , x2 := 1.6 ]
reg_main_out[ asset_type == "Illiquid assets" & grp == "Assets low | Income low" , x2 := 2.1 ]

reg_main_out[ asset_type == "Liquid assets" , asset_type2 := "Liquid" ]
reg_main_out[ asset_type == "Illiquid assets" , asset_type2 := "Illiquid" ]
reg_main_out[ , asset_type2 := factor(asset_type2, levels=c("Liquid","Illiquid")) ]

cairo_pdf("output/fig_4.pdf", width=7, height=2.75)
(plot <- (ggplot(data=reg_main_out, aes(y=coef, x=x2, colour=asset_type2, shape=asset_type2, fill=asset_type2))
          + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", size=.5)
          + geom_errorbar(aes(ymin=lower95, ymax=upper95), width=0, size=0.4, show.legend = F)  
          + geom_errorbar(aes(ymin=lower90, ymax=upper90), width=0, size=0.8)  
          + geom_point(size=2.8)
          + facet_wrap( ~ dv2)
          + coord_flip()
          + ylab("Marginal effect")
          + theme_bw()
          + theme(legend.position="right", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                  panel.border = element_rect(colour="gray"))
          + guides(colour = guide_legend(ncol = 1))
          + scale_x_continuous(labels=c("Income-buffered\n(income high/\nassets low)",
                                        "Asset-buffered\n(income low/\nassets high)",
                                        "Economically precarious\n(income low/\nassets low)"),
                               breaks=c(1.05,1.55,2.05), limits=c(0.8,2.3))
          + ylim(-0.02,0.151)
          + scale_color_manual(name="Asset type", values=c("darkgreen","black"))
          + scale_fill_manual(name="Asset type", values=c("darkgreen","black"))
          + scale_shape_manual(name="Asset type", values=c(21,22))
))
dev.off()




#########################
### Interaction model ###
#########################

mod2.1 <- lm(gvt_help_unemp_num_norm ~ tot_net_inc_USD_v1_dec_v2*liquid_assets_USD_v1_dec + age + age_sqr + gender + highest_educ_v2 + unemployed +
                student + retired + employed + homeowner + no_chld + married + country_name, cs_survey, weights=cs_survey$wgt)

mod2.2 <- lm(gvt_help_unemp_num_norm ~ tot_net_inc_USD_v1_dec_v2*ass_real_estate_USD_v1_dec + age + age_sqr + gender + highest_educ_v2 + unemployed +
                student + retired + employed + homeowner + no_chld + married + country_name, cs_survey, weights=cs_survey$wgt)

mod2.3 <- lm(gvt_reduce_inc_diff_num_norm ~ tot_net_inc_USD_v1_dec_v2*liquid_assets_USD_v1_dec + age + age_sqr + gender + highest_educ_v2 + unemployed +
                student + retired + employed + homeowner + no_chld + married + country_name, cs_survey, weights=cs_survey$wgt)

mod2.4 <- lm(gvt_reduce_inc_diff_num_norm ~ tot_net_inc_USD_v1_dec_v2*ass_real_estate_USD_v1_dec + age + age_sqr + gender + highest_educ_v2 + unemployed +
                student + retired + employed + homeowner + no_chld + married + country_name, cs_survey, weights=cs_survey$wgt)

### Appendix Table E.1
lab_covar2 <- c("Income","Liquid assets","Illiquid  assets",
                "Age","Age square","Male","Education: high school","Education: vocational training",
                "Education: college", "Education: postgrad",
                "Unemployed","Student","Retired","Employed","Homeowner","Number of children","Married",
                "Income $\\times$ Liquid assets","Income $\\times$ Illiquid assets")

stargazer(mod2.1, mod2.2, mod2.3, mod2.4,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3,
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar2,
          add.lines = list(c("Mean DV",
                             round(mean(mod2.1$fitted.values),3),
                             round(mean(mod2.2$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$")))

## Appendix Figure E.1
d_plot_mod2.1 <- cbind(as.data.table(ggplot_build(plot_model(mod2.1, type="int"))$data[[1]]),
                        as.data.table(ggplot_build(plot_model(mod2.1, type="int"))$data[[2]])[,.(ymax,ymin) ])

d_plot_mod2.2 <- cbind(as.data.table(ggplot_build(plot_model(mod2.2, type="int"))$data[[1]]),
                        as.data.table(ggplot_build(plot_model(mod2.2, type="int"))$data[[2]])[,.(ymax,ymin) ])

d_plot_mod2.3 <- cbind(as.data.table(ggplot_build(plot_model(mod2.3, type="int"))$data[[1]]),
                       as.data.table(ggplot_build(plot_model(mod2.3, type="int"))$data[[2]])[,.(ymax,ymin) ])

d_plot_mod2.4 <- cbind(as.data.table(ggplot_build(plot_model(mod2.4, type="int"))$data[[1]]),
                       as.data.table(ggplot_build(plot_model(mod2.4, type="int"))$data[[2]])[,.(ymax,ymin) ])

d_plot_mod2.1[ , type := "Liquid assets" ]
d_plot_mod2.2[ , type := "Illiquid assets" ]

d_plot1 <- rbind(d_plot_mod2.1, d_plot_mod2.2)
d_plot1[ ,  dv_label := "Support for UI insurance" ]

d_plot_mod2.3[ , type := "Liquid assets" ]
d_plot_mod2.4[ ,type := "Illiquid assets" ]

d_plot2 <- rbind(d_plot_mod2.3, d_plot_mod2.4)
d_plot2[ , dv_label := "Support for redistribution" ]

d_plot_int_comb <- as.data.table(rbind(d_plot1, d_plot2))

d_plot_int_comb[ group == 1, ass_grp := "bottom" ]
d_plot_int_comb[ group == 2, ass_grp := "top" ]

d_plot_int_comb[ , dv_label := factor(dv_label, levels=c("Support for UI insurance","Support for redistribution")) ]
d_plot_int_comb[ , type := factor(type, levels=c("Liquid assets","Illiquid assets"))]

(p1.1 <- (ggplot(data=d_plot_int_comb[ dv_label == "Support for UI insurance" ],
                 aes(y=y, x=x, colour=ass_grp, linetype=ass_grp))
          + geom_ribbon(aes(ymin=ymin, ymax=ymax, fill=ass_grp), alpha=0.25, colour=NA)
          + geom_line(size=1)
          + facet_wrap(~ type)
          + theme_bw()
          + xlab("Net income (in deciles)")
          + ylab("Predicted support")
          + theme(legend.position="right", panel.grid.minor = element_blank(),
                  panel.border = element_rect(colour="gray"))
          + guides(colour = guide_legend(ncol = 1))
          + coord_cartesian(ylim=c(0.52,0.85))
          + scale_y_continuous(breaks=seq(0.5,0.8,0.1))
          + scale_x_continuous(breaks=seq(2,10,2))
          + scale_color_manual(name="Asset group", values=c("black","darkgreen"))
          + scale_fill_manual(name="Asset group", values=c("black","darkgreen"))
          + scale_linetype_manual(name="Asset group", values=c("dotdash","solid"))
          + ggtitle("Support for unemployment insurance")
))

(p1.2 <- (ggplot(data=d_plot_int_comb[ dv_label == "Support for redistribution" ],
                 aes(y=y, x=x, colour=ass_grp, linetype=ass_grp))
          + geom_ribbon(aes(ymin=ymin, ymax=ymax, fill=ass_grp), alpha=0.25, colour=NA)
          + geom_line(size=1)
          + facet_wrap(~ type)
          + theme_bw()
          + xlab("Net income (in deciles)")
          + ylab("Predicted support")
          + theme(legend.position="right", panel.grid.minor = element_blank(),
                  panel.border = element_rect(colour="gray"))
          + guides(colour = guide_legend(ncol = 1))
          + coord_cartesian(ylim=c(0.52,0.85))
          + scale_y_continuous(breaks=seq(0.5,0.8,0.1))
          + scale_x_continuous(breaks=seq(2,10,2))
          + scale_color_manual(name="Asset group", values=c("black","darkgreen"))
          + scale_fill_manual(name="Asset group", values=c("black","darkgreen"))
          + scale_linetype_manual(name="Asset group", values=c("dotdash","solid"))
          + ggtitle("Support for redistribution")
))


cairo_pdf("output/si_fig_e1.pdf",width=6, height=6.5)
grid.arrange(p1.1,p1.2)
dev.off()



#########################
### Robustness checks ###
#########################

### Remove homeowner indicator
mod1.1_R1 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.2_R1 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + no_chld + married
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.3_R1 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.4_R1 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + no_chld + married 
               | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

### Appendix Table D.2
lab_covar_R1 <- c("Asset-buffered","Income-buffered","Economically-precarious",
                "Asset-buffered","Income-buffered","Economically-precarious",
                "Age","Age square","Male","Education: high school","Education: vocational training",
                "Education: college", "Education: postgrad",
                "Unemployed","Student","Retired","Employed","Number of children","Married")

stargazer(mod1.1_R1,mod1.3_R1,mod1.2_R1,mod1.4_R1,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3,  
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar_R1,
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1_R1$fitted.values),3), 
                             round(mean(mod1.3_R1$fitted.values),3), 
                             round(mean(mod1.2_R1$fitted.values),3), 
                             round(mean(mod1.4_R1$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$")))




### Subset to those how are certain about the assets
cs_survey_certain <- cs_survey[ how_certain_asset_q %in% c(2,3) ]

mod1.1_R2 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married | country_name | 0 | country_name, cs_survey_certain, weights = cs_survey_certain$wgt )

mod1.2_R2 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married | country_name | 0 | country_name, cs_survey_certain, weights = cs_survey_certain$wgt )

mod1.3_R2 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married | country_name | 0 | country_name, cs_survey_certain, weights = cs_survey_certain$wgt )

mod1.4_R2 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married | country_name | 0 | country_name, cs_survey_certain, weights = cs_survey_certain$wgt )

### Appendix Table D.4 
lab_covar_R2 <- c("Asset-buffered","Income-buffered","Economically-precarious",
                   "Asset-buffered","Income-buffered","Economically-precarious",
                   "Age","Age square","Male","Education: high school","Education: vocational training",
                   "Education: college", "Education: postgrad",
                   "Unemployed","Student","Retired","Employed","Homeowner","Number of children","Married")

stargazer(mod1.1_R2, mod1.3_R2, mod1.2_R2, mod1.4_R2,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3, 
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar_R2, 
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1_R2$fitted.values),3), 
                             round(mean(mod1.2_R2$fitted.values),3), 
                             round(mean(mod1.3_R2$fitted.values),3), 
                             round(mean(mod1.4_R2$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$")))

# Results
out2 <- data.table(rbind(cbind(mod1.1_R2$coefficients[1:3], confint(mod1.1_R2, level=0.95)[1:3,], confint(mod1.1_R2, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Liquid assets", grp=rownames(mod1.1_R2$coefficients)[1:3]),
                         cbind(mod1.2_R2$coefficients[1:3], confint(mod1.2_R2, level=0.95)[1:3,], confint(mod1.2_R2, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Liquid assets", grp=rownames(mod1.2_R2$coefficients)[1:3]),
                         cbind(mod1.3_R2$coefficients[1:3], confint(mod1.3_R2, level=0.95)[1:3,], confint(mod1.3_R2, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.3_R2$coefficients)[1:3]),
                         cbind(mod1.4_R2$coefficients[1:3], confint(mod1.4_R2, level=0.95)[1:3,], confint(mod1.4_R2, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.4_R2$coefficients)[1:3])))

out2[ , 1:5 := lapply( .SD, function (x) as.numeric(as.character(x)) ), .SDcols=c(1:5) ]

out2[ , covars := "(3) Subset: high certainty"]



### With controls for income loss, risk attitudes, and economic worries
mod1.1_R3 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + take_risk + inc_expect_num + worried_overall_econ_sec | country_name | 0 | country_name, cs_survey,
                   weights = cs_survey$wgt)

mod1.2_R3 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + take_risk + inc_expect_num + worried_overall_econ_sec | country_name | 0 | country_name, cs_survey,
                   weights = cs_survey$wgt)

mod1.3_R3 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + take_risk + inc_expect_num + worried_overall_econ_sec | country_name | 0 | country_name, cs_survey,
                   weights = cs_survey$wgt)

mod1.4_R3 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + take_risk + inc_expect_num + worried_overall_econ_sec | country_name | 0 | country_name, cs_survey,
                   weights = cs_survey$wgt)


### Appendix Table D.3 
lab_covar_R3 <- c("Asset-buffered","Income-buffered","Economically-precarious",
                   "Asset-buffered","Income-buffered","Economically-precarious",
                   "Age","Age square","Male","Education: high school","Education: vocational training",
                   "Education: college", "Education: postgrad",
                   "Unemployed","Student","Retired","Employed","Homeowner","Number of children","Married",
                   "Risk aversion","Income expectations","Worried economic security")

stargazer(mod1.1_R3, mod1.3_R3, mod1.2_R3, mod1.4_R3,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3, label="tab:mod_UI_rob1x", 
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar_R3, 
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1_R3$fitted.values),3), 
                             round(mean(mod1.2_R3$fitted.values),3), 
                             round(mean(mod1.3_R3$fitted.values),3), 
                             round(mean(mod1.4_R3$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$")))

# Results
out3 <- data.table(rbind(cbind(mod1.1_R3$coefficients[1:3], confint(mod1.1_R3, level=0.95)[1:3,], confint(mod1.1_R3, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Liquid assets", grp=rownames(mod1.1_R3$coefficients)[1:3]),
                         cbind(mod1.2_R3$coefficients[1:3], confint(mod1.2_R3, level=0.95)[1:3,], confint(mod1.2_R3, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Liquid assets", grp=rownames(mod1.2_R3$coefficients)[1:3]),
                         cbind(mod1.3_R3$coefficients[1:3], confint(mod1.3_R3, level=0.95)[1:3,], confint(mod1.3_R3, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.3_R3$coefficients)[1:3]),
                         cbind(mod1.4_R3$coefficients[1:3], confint(mod1.4_R3, level=0.95)[1:3,], confint(mod1.4_R3, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.4_R3$coefficients)[1:3])))

out3[ , 1:5 := lapply( .SD, function (x) as.numeric(as.character(x)) ), .SDcols=c(1:5) ]

out3[ , covars := "(2) + risk covariates"]


### With country covariates 
load("data/cpds_2017.RData")

cs_survey[ , country := country_name ]
cs_survey[ country_name == "England", country := "United Kingdom"]

cs_survey_cty <- merge(cs_survey, cpds_2017, by="country", all.x=T)

# Models 
mod1.1_R4 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + unemp + realgdpgr + sstran + elderly + gov_right3 + prop + fed, cs_survey_cty, weights = cs_survey_cty$wgt)

mod1.2_R4 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + unemp + realgdpgr + sstran + elderly + gov_right3 + prop + fed, cs_survey_cty, weights = cs_survey_cty$wgt)

mod1.3_R4 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + unemp + realgdpgr + sstran + elderly + gov_right3 + prop + fed, cs_survey_cty, weights = cs_survey_cty$wgt)

mod1.4_R4 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired +  + 
                     employed + homeowner + no_chld + married + unemp + realgdpgr + sstran + elderly + gov_right3 + prop + fed, cs_survey_cty, weights = cs_survey_cty$wgt)


### Appendix Table D.5
lab_covar_R4 <- c("Asset-buffered","Income-buffered","Economically-precarious",
                   "Asset-buffered","Income-buffered","Economically-precarious",
                   "Age","Age square","Male","Education: high school","Education: vocational training",
                   "Education: college", "Education: postgrad",
                   "Unemployed","Student","Retired","Employed","Homeowner","Number of children","Married",
                   "Unemployment rate","Real GDP growth","Social spending","Share elderly",
                   "Seat share right parties","Modified prop. representation","Prop. representation","Strong federalism")

stargazer(mod1.1_R4, mod1.3_R4, mod1.2_R4, mod1.4_R4,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3, label="tab:mod_UI_rob1x", 
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
          covariate.labels = lab_covar_R4, 
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1_R4$fitted.values),3), 
                             round(mean(mod1.2_R4$fitted.values),3), 
                             round(mean(mod1.3_R4$fitted.values),3), 
                             round(mean(mod1.4_R4$fitted.values),3))))

# Results
out4 <- data.table(rbind(cbind(mod1.1_R4$coefficients[2:4], confint(mod1.1_R4, level=0.95)[2:4,], confint(mod1.1_R4, level=0.9)[2:4,], dv="gvt_help_unemp_num_norm", asset_type="Liquid assets", grp=rownames(mod1.1_R4$coefficients)[2:4]),
                         cbind(mod1.2_R4$coefficients[2:4], confint(mod1.2_R4, level=0.95)[2:4,], confint(mod1.2_R4, level=0.9)[2:4,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Liquid assets", grp=rownames(mod1.2_R4$coefficients)[2:4]),
                         cbind(mod1.3_R4$coefficients[2:4], confint(mod1.3_R4, level=0.95)[2:4,], confint(mod1.3_R4, level=0.9)[2:4,], dv="gvt_help_unemp_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.3_R4$coefficients)[2:4]),
                         cbind(mod1.4_R4$coefficients[2:4], confint(mod1.4_R4, level=0.95)[2:4,], confint(mod1.4_R4, level=0.9)[2:4,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.4_R4$coefficients)[2:4])))

out4[ , 1:5 := lapply( .SD, function (x) as.numeric(as.character(x)) ), .SDcols=c(1:5) ]

out4[ , covars := "(4) + country covariates"]



### With probability of unemployment
cs_survey[ , prob_unemp_pct := prob_unemp / 100  ]

mod1.1_R5 <- felm(gvt_help_unemp_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + 
                 no_chld + married + prob_unemp_pct | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.2_R5 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_liq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + 
                 married + prob_unemp_pct | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.3_R5 <- felm(gvt_help_unemp_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + 
                 married + prob_unemp_pct | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

mod1.4_R5 <- felm(gvt_reduce_inc_diff_num_norm ~ grp_med_illiq + age + age_sqr + gender + highest_educ_v2 + unemployed + student + retired + employed + homeowner + no_chld + 
                 married + prob_unemp_pct | country_name | 0 | country_name, cs_survey, weights = cs_survey$wgt)

### Appendix Table D.6
lab_covar_R5 <- c("Asset-buffered","Income-buffered","Economically-precarious",
                "Asset-buffered","Income-buffered","Economically-precarious",
                "Age","Age square","Male","Education: high school","Education: vocational training",
                "Education: college", "Education: postgradudate",
                "Unemployed","Student","Retired","Employed","Homeowner","Number of children",
                "Married","Probability of unemployment")

stargazer(mod1.1_R5,mod1.3_R5,mod1.2_R5,mod1.4_R5,
          omit=c("factor","Constant","country_name"), digits.extra=0, digits=3,  
          align=T, no.space=T, omit.stat=c("ser","F"), type="text", column.sep.width = "-10pt",
          dep.var.labels = c("Unemployment Insurance","Redistribution"),
            covariate.labels = lab_covar_R5,
          add.lines = list(c("Mean DV",
                             round(mean(mod1.1_R5$fitted.values),3), 
                             round(mean(mod1.3_R5$fitted.values),3), 
                             round(mean(mod1.2_R5$fitted.values),3), 
                             round(mean(mod1.4_R5$fitted.values),3)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$")))




### Figure 5
reg_out_comb <- data.table(rbind(cbind(mod1.1$coefficients[1:3], confint(mod1.1, level=0.95)[1:3,], confint(mod1.1, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Liquid assets", grp=rownames(mod1.1$coefficients)[1:3]),
                             cbind(mod1.2$coefficients[1:3], confint(mod1.2, level=0.95)[1:3,], confint(mod1.2, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Liquid assets", grp=rownames(mod1.2$coefficients)[1:3]),
                             cbind(mod1.3$coefficients[1:3], confint(mod1.3, level=0.95)[1:3,], confint(mod1.3, level=0.9)[1:3,], dv="gvt_help_unemp_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.3$coefficients)[1:3]),
                             cbind(mod1.4$coefficients[1:3], confint(mod1.4, level=0.95)[1:3,], confint(mod1.4, level=0.9)[1:3,], dv="gvt_reduce_inc_diff_num_norm", asset_type="Illiquid assets", grp=rownames(mod1.4$coefficients)[1:3])))

reg_out_comb[ , 1:5 := lapply( .SD, function (x) as.numeric(as.character(x)) ), .SDcols=c(1:5) ]

reg_out_comb[ , covars := "(1) Baseline"]

reg_rob_out_all <- rbind(reg_out_comb,out2,out3)

names(reg_rob_out_all) <- c("coef","lower95","upper95","lower90","upper90","dv","asset_type","grp","covars")

reg_rob_out_all[ , grp := gsub("grp_med_liqLiq ","",grp)]
reg_rob_out_all[ , grp := gsub("grp_med_illiqIlliq ","",grp)]

reg_rob_out_all[ , asset_type := factor(asset_type, levels=c("Liquid assets","Illiquid assets")) ]

reg_rob_out_all[ dv == "gvt_help_unemp_num_norm" , dv2 := "Support for unemployment insurance" ]
reg_rob_out_all[ dv == "gvt_reduce_inc_diff_num_norm" , dv2 := "Support for redistribution" ]

reg_rob_out_all[ , covars := factor(covars, levels=c("(1) Baseline","(2) + risk covariates","(3) Subset: high certainty","(4) + country covariates")) ]

reg_rob_out_all[ grp == "Assets low | Income high" & covars == "(3) Subset: high certainty", x3 := 1 ]
reg_rob_out_all[ grp == "Assets low | Income high" & covars == "(2) + risk covariates", x3 := 1.1 ]
reg_rob_out_all[ grp == "Assets low | Income high" & covars == "(1) Baseline", x3 := 1.2 ]

reg_rob_out_all[ grp == "Assets high | Income low" & covars == "(3) Subset: high certainty", x3 := 1.7 ]
reg_rob_out_all[ grp == "Assets high | Income low" & covars == "(2) + risk covariates", x3 := 1.8 ]
reg_rob_out_all[ grp == "Assets high | Income low" & covars == "(1) Baseline", x3 := 1.9 ]

reg_rob_out_all[ grp == "Assets low | Income low" & covars == "(3) Subset: high certainty", x3 := 2.4 ]
reg_rob_out_all[ grp == "Assets low | Income low" & covars == "(2) + risk covariates", x3 := 2.5 ]
reg_rob_out_all[ grp == "Assets low | Income low" & covars == "(1) Baseline", x3 := 2.6 ]

(plot2.1 <- (ggplot(data=reg_rob_out_all[ dv2 == "Support for unemployment insurance" & covars != "(4) + country covariates" ], aes(y=coef, x=x3, fill=covars, shape=covars, colour=covars))
             + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", size=.5)
             + geom_errorbar(aes(ymin=lower95, ymax=upper95), width=0, size=0.4, show.legend = F)  
             + geom_errorbar(aes(ymin=lower90, ymax=upper90), width=0, size=0.8)  
             + geom_point(size=2.1)
             + facet_grid( ~ asset_type)
             + coord_flip()
             + ylab("Marginal effect")
             + theme_bw()
             + theme(legend.position="right", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                     panel.border = element_rect(colour="gray"))
             + guides(colour = guide_legend(ncol = 1))
             + scale_x_continuous(labels=c("Asset-buffered\n(income low/\nassets high)",
                                           "Income-buffered\n(income high/\nassets low)",
                                           "Economically precarious\n(income low/\nassets low)"),
                                  breaks=c(1.1,1.8,2.5), limits=c(0.9,2.7))
             + ylim(-0.02,0.151)
             + scale_color_manual(name="Model", values=c("black","darkgreen","blue","red"))
             + scale_fill_manual(name="Model", values=c("black","darkgreen","blue","red"))
             + scale_shape_manual(name="Model", values=c(21:24))
             + ggtitle("Support for unemployment insurance")
))

(plot2.2 <- (ggplot(data=reg_rob_out_all[ dv2 == "Support for redistribution" & covars != "(4) + country covariates" ], aes(y=coef, x=x3, fill=covars, shape=covars, colour=covars))
             + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", size=.5)
             + geom_errorbar(aes(ymin=lower95, ymax=upper95), width=0, size=0.4, show.legend = F)  
             + geom_errorbar(aes(ymin=lower90, ymax=upper90), width=0, size=0.8)  
             + geom_point(size=2.1)
             + facet_grid( ~ asset_type)
             + coord_flip()
             + ylab("Marginal effect")
             + theme_bw()
             + theme(legend.position="right", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                     panel.border = element_rect(colour="gray"))
             + guides(colour = guide_legend(ncol = 1))
             + scale_x_continuous(labels=c("Income-buffered\n(income high/\nassets low)",
                                           "Asset-buffered\n(income low/\nassets high)",
                                           "Economically precarious\n(income low/\nassets low)"),
                                  breaks=c(1.1,1.8,2.5), limits=c(0.9,2.7))
             + ylim(-0.02,0.151)
             + scale_color_manual(name="Model", values=c("black","darkgreen","blue","red"))
             + scale_fill_manual(name="Model", values=c("black","darkgreen","blue","red"))
             + scale_shape_manual(name="Model", values=c(21:24))
             + ggtitle("Support for redistribution")
))

cairo_pdf("output/fig_5.pdf", width=8, height=6)
ggarrange(plot2.1, NULL, plot2.2, ncol=1, nrow=3, heights = c(1,0.05,1), common.legend = TRUE, legend="right")
dev.off()




################################################
### Correlation figures with decile measures ###
################################################

# Calculate correlation for liquid wealth
inc_liq <- na.omit(cs_survey[ ,.(country_name,tot_net_inc_USD_v1_dec_v2,liquid_assets_USD_v1_dec,wgt)])

corr_cty <- inc_liq[ , weightedCorr(tot_net_inc_USD_v1_dec_v2, liquid_assets_USD_v1_dec, method=c("Pearson"), weights = wgt), by=c("country_name") ]

names(corr_cty)[2] <- "corr_inc_liq_ass"

# Illiquid wealth
inc_illiq <- na.omit(cs_survey[ ,.(country_name,tot_net_inc_USD_v1_dec_v2,ass_real_estate_USD_v1_dec,wgt)])

corr_cty_illiquid <- inc_illiq[ , weightedCorr(tot_net_inc_USD_v1_dec_v2, ass_real_estate_USD_v1_dec, method=c("Pearson"), weights = wgt), by=c("country_name") ]

names(corr_cty_illiquid)[2] <- "corr_inc_illiq_ass"


### Figure 7 
cairo_pdf("output/fig_7.pdf", width=4.5, height=3)
(plot <- (ggplot(data=corr_cty, aes(y=corr_inc_liq_ass, x=reorder(country_name,corr_inc_liq_ass)))
          + geom_bar(stat = "identity", position = "dodge", width=0.7, fill="darkgreen")
          + geom_text(aes(label=round(corr_inc_liq_ass,2)), size=2.5, hjust=1.25, colour="white")
          + theme_bw()
          + ylab("Correlation(income, liquid assets)")
          + coord_flip()
          + scale_y_continuous(breaks = seq(0,0.6,0.2), limits = c(0,0.6))
          + theme(legend.position="bottom", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                  panel.border = element_rect(colour="gray"), panel.grid.major.x = element_blank())
))
dev.off()



### Binscatter plot by country 
cty_list <- sort(unique(cs_survey$country_name))

bins_out_liq_list <- list()
bins_out_illiq_list <- list()
fit.liq.d.list <- list()
fit.illiq.d.list <- list()

for(i in 1:9){
  d_cty <- cs_survey[ country_name == cty_list[i] ]
  
  # Liquid
  bins_plot_liq <- binsreg(d_cty$liquid_assets_USD_v1, d_cty$tot_net_inc_USD_v1, weights = d_cty$wgt)
  
  bins_out_liq <- as.data.table(bins_plot_liq$data.plot$`Group Full Sample`$data.dots)
  
  bins_out_liq[ , cty := cty_list[i] ]
  
  bins_out_liq_list[[i]] <- bins_out_liq
  
  # Illiquid
  bins_plot_illiq <- binsreg(d_cty$ass_real_estate_USD_v1, d_cty$tot_net_inc_USD_v1, weights = d_cty$wgt)
  
  bins_out_illiq <- as.data.table(bins_plot_illiq$data.plot$`Group Full Sample`$data.dots)
  
  bins_out_illiq[ , cty := cty_list[i] ]
  
  bins_out_illiq_list[[i]] <- bins_out_illiq
  
  # Get lm output
  fit.liq <- lm(fit ~ x, bins_out_liq)
  
  fit.liq.d <- as.data.table(fit.liq$coefficients[2])
  fit.liq.d[ , cty := cty_list[i] ]
  fit.liq.d.list[[i]] <- fit.liq.d
  
  fit.illiq <- lm(fit ~ x, bins_out_illiq)
  
  fit.illiq.d <- as.data.table(fit.illiq$coefficients[2])
  fit.illiq.d[ , cty := cty_list[i] ]
  fit.illiq.d.list[[i]] <- fit.illiq.d
  
  print(i)
}

bins_out_illiq_comb <- rbindlist(bins_out_illiq_list)
bins_out_liq_comb <- rbindlist(bins_out_liq_list)

fit_illiq_comb <- rbindlist(fit.illiq.d.list)
fit_liq_comb <- rbindlist(fit.liq.d.list)

bins_out_illiq_comb[ , type := "illiquid"]
bins_out_liq_comb[ , type := "liquid"]

bins_out_comb <- rbind(bins_out_illiq_comb, bins_out_liq_comb)

fit_illiq_comb[ , type := "illiquid"]
fit_liq_comb[ , type := "liquid"]

fit_out_comb <- rbind(fit_illiq_comb, fit_liq_comb)

fit_out_comb[ type == "liquid" , fit := paste("Liquid: \u03b2=",round(V1,2), sep="")]
fit_out_comb[ type == "illiquid" , fit := paste("Illiquid: \u03b2=",round(V1,2), sep="")]

fit_out_comb[ type == "liquid" , fit2 := paste("\u03b2=",round(V1,2), sep="")]
fit_out_comb[ type == "illiquid" , fit2 := paste("\u03b2=",round(V1,2), sep="")]

bins_out_comb[ , fit_tsd := fit/1000 ]
bins_out_comb[ , x_tsd := x/1000 ]

bins_out_comb[ , cty := factor(cty,levels=sort(as.character(cty_list)))]


### Figure 6 
cairo_pdf("output/fig_6.pdf", width=7, height=5)
(plot <- (ggplot(data=bins_out_comb, aes(y=fit_tsd, x=x_tsd, fill=type, color=type, shape=type))
          + geom_point(alpha=0.5)
          + stat_smooth(method="lm")
          + geom_text(data=fit_out_comb[ type == "illiquid" ], mapping = aes(x = 210, y = 50, label = fit2), hjust = 1, vjust = 0, size=2.5)
          + geom_text(data=fit_out_comb[ type == "liquid" ], mapping = aes(x = 210, y = 0, label = fit2), hjust = 1, vjust = 0, size=2.5)
          + theme_bw()
          + ylab("Assets (thsd. USD)")
          + xlab("Income (thsd. USD)")
          + facet_wrap(~ cty, ncol=5)
          + scale_color_manual(name="Asset type", values=c("black","darkgreen"))
          + scale_fill_manual(name="Asset type", values=c("black","darkgreen"))
          + scale_shape_manual(name="Asset type", values=c(21,22))
          + theme(legend.position="bottom", panel.grid.minor = element_blank(), 
                  panel.border = element_rect(colour="gray"), panel.grid.major.x = element_blank())
))
dev.off()




### Figure 8 
d_svy_wgt2 <- cs_survey[ !is.na(grp_med_liq) ] %>% as_survey_design(weights = wgt)

d_grp_wgt <- as.data.table(d_svy_wgt2 %>% group_by(country_name, grp_med_liq) %>% summarise(pct = survey_prop()))

d_grp_wgt[ grp_med_liq == "Liq Assets high | Income low", type := "Asset-buffered:\n(low income &\nhigh liquid assets)" ]
d_grp_wgt[ grp_med_liq == "Liq Assets low | Income high", type := "Income-buffered:\n(high income &\nlow liquid assets)" ]
d_grp_wgt[ grp_med_liq == "Liq Assets high | Income high", type := "Truly wealthy:\n(high income &\nhigh liquid assets)" ]
d_grp_wgt[ grp_med_liq == "Liq Assets low | Income low", type := "Economically precarious:\n(low income &\nlow liquid assets)" ]

d_grp_wgt.1 <- d_grp_wgt[ grp_med_liq %in% c("Liq Assets high | Income low","Liq Assets low | Income high") ] 

d_grp_wgt.1[ , country_name := factor(country_name, levels=c("USA","Spain","Germany","Sweden","Canada","Netherlands","France","England","Denmark")) ]
d_grp_wgt.1[ , pct := pct*100 ]

cairo_pdf("output/fig_8.pdf", width=6, height=3)
(plot <- (ggplot(data=d_grp_wgt.1, aes(y=pct, x=country_name, fill=type, label=round(pct,1)))
          + geom_bar(stat="identity")
          + geom_text(colour = "white", size = 2.5, position = position_stack(vjust = 0.5))
          + theme_bw()
          + ylab("Population shares (%)")
          + xlab("Income group")
          + coord_flip()
          + theme(legend.position="right", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                  panel.border = element_rect(colour="gray"), plot.title = element_text(hjust = 0.5))
          + guides(fill = guide_legend(nrow = 2, keyheight=2.5))
          + scale_fill_manual("Group",values=c("blue","black"))
          + ylim(0,40)
))
dev.off()




########################################
### Correlation with social spending ###
########################################

### OECD SOCX
load("data/oecd_socx.RData") # socx

### Figure 9

# Subset without health and pensions
socx_v1 <- socx[ UNIT == "PCT_GDP" & SOURCE == "10" ]

# Sum
socx_v1 <- socx_v1[ , lapply(.SD, function(x) sum(x, na.rm=T) ), by=c("country_name"), .SDcols= "Value" ]

corr_cty_spend <- merge(corr_cty, socx_v1, by=c("country_name"))
corr_cty_spend <- merge(corr_cty_spend, corr_cty_illiquid, by=c("country_name"))

cor_test <- cor.test(corr_cty_spend$Value, corr_cty_spend$corr_inc_liq_ass, method=c("pearson"))

cor_out <- paste("\u03b2 =",round(cor_test$estimate,2),", p =",round(cor_test$p.value,3))
cor_out <- paste("r = ",round(cor_test$estimate,2),", p = ",round(cor_test$p.value,3),sep="")

cairo_pdf("output/fig_9.pdf", width=5.5, height=4)
(plot <- (ggplot(data=corr_cty_spend, aes(y=Value, x=corr_inc_liq_ass))
          + stat_smooth(method="lm", size=1.25)
          + geom_point(size=2)
          + geom_text_repel(aes(label=country_name), size=3)
          + annotate("text", x = 0.32, y = 1, label = cor_out, hjust = 0.5, vjust = 0, size=2.25)
          + theme_bw()
          + ylab("Government social spending (% of GDP)")
          + xlab("Correlation(income,liquid assets)")
          + theme(legend.position="bottom", panel.grid.minor = element_blank(), 
                  panel.border = element_rect(colour="gray"), panel.grid.major.x = element_blank())
          + ylim(1,16)
))
dev.off()



### Appendix Figure G.1

# Share of government expenditures
socx_v2 <- socx[ UNIT == "PCT_GOV" ]

socx_v2 <- socx_v2[ , lapply(.SD, function(x) sum(x, na.rm=T) ), by=c("country_name"), .SDcols= "Value" ]

corr_cty_spend_share <- merge(corr_cty, socx_v2, by=c("country_name"))

cor_test2 <- cor.test(corr_cty_spend_share$Value, corr_cty_spend_share$corr_inc_liq_ass, method=c("pearson"))

cor_out2 <- paste("\u03b2 =",round(cor_test2$estimate,2),", p =",round(cor_test2$p.value,3))
cor_out2 <- paste("r = ",round(cor_test2$estimate,2),", p = ",round(cor_test2$p.value,3),sep="")

cairo_pdf("output/si_fig_g1.pdf", width=5.5, height=4)
(plot <- (ggplot(data=corr_cty_spend_share, aes(y=Value, x=corr_inc_liq_ass))
          + stat_smooth(method="lm", size=1.25)
          + geom_point(size=2)
          + geom_text_repel(aes(label=country_name), size=3)
          + annotate("text", x = 0.32, y = 7, label = cor_out2, hjust = 0.5, vjust = 0, size=2.25)
          + theme_bw()
          + ylab("Government social spending (% of total spending)")
          + xlab("Correlation(income,liquid assets)")
          + theme(legend.position="bottom", panel.grid.minor = element_blank(), 
                  panel.border = element_rect(colour="gray"), panel.grid.major.x = element_blank())
))
dev.off()



### Appendix Figure F.1
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(1,2,3) , inc_grp_v3 := "income: low" ]
cs_survey[ tot_net_inc_USD_v1_dec_v2 %in% c(8,9,10) , inc_grp_v3 := "income: high" ]
cs_survey[ liquid_assets_USD_v1_dec %in% c(1,2,3) , ass_grp_v3 := "liquid assets: low" ]
cs_survey[ liquid_assets_USD_v1_dec %in% c(8,9,10) , ass_grp_v3 := "liquid assets: high" ]

# Off-Diagonals
cs_survey[ ass_grp_v3 == "liquid assets: high" & inc_grp_v3 == "income: low", type := "Asset-buffered:\n(low income &\nhigh liquid assets)" ]
cs_survey[ ass_grp_v3 == "liquid assets: low" & inc_grp_v3 == "income: high", type := "Income-buffered:\n(high income &\nlow liquid assets)" ]

d_svy_wgt3 <- cs_survey[ !is.na(inc_grp_v3) & !is.na(ass_grp_v3) ] %>% as_survey_design(weights = wgt)

grp_share <- as.data.table(d_svy_wgt3 %>% group_by(country_name, type) %>% summarise(pct = survey_prop()))
grp_share[ , pct := pct * 100 ]

grp_share[ , country_name := factor(country_name, levels=c("USA","Spain","Canada","Netherlands","France","Sweden","Germany","England","Denmark")) ]

cairo_pdf("output/si_fig_f1.pdf", width=6, height=3)
(plot <- (ggplot(data=grp_share[ !is.na(type) ], aes(y=pct, x=country_name, fill=type, label=round(pct,1)))
          + geom_bar(stat="identity")
          + geom_text(colour = "white", size = 2.5, position = position_stack(vjust = 0.5))
          + theme_bw()
          + ylab("Population shares (%)")
          + xlab("Income group")
          + coord_flip()
          + theme(legend.position="right", panel.grid.minor = element_blank(), axis.title.y = element_blank(),
                  panel.border = element_rect(colour="gray"), plot.title = element_text(hjust = 0.5))
          + guides(fill = guide_legend(nrow = 2, keyheight=2.5))
          + scale_fill_manual("Group",values=c("blue","black"))
))
dev.off()






